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ABSTRACT 

We performed a series of high-resolution (up to 1024'^) direct numerical simulations of hydro and 
MHD turbulence. Our simulations correspond to the "strong" MHD turbulence regime that cannot 
be treated perturbatively. We found that for simulations with normal viscosity the slopes for en- 
ergy spectra of MHD are similar to ones in hydro, although slightly more shallower. However, for 
simulations with hyper viscosity the slopes were very different, for instance, the slopes for hydro sim- 
ulations showed pronounced and well-defined bottleneck effect, while the MHD slopes were relatively 
much less affected. We believe that this is indicative of MHD strong turbulence being less local than 
Kolmogorov turbulence. This calls for revision of MHD strong turbulence models that assume local 
"as-in-hydro case" cascading. Non-locality of MHD turbulence casts doubt on numerical determina- 
tion of the slopes with currently available (512^-1024'^) numerical resolutions, including simulations 
with normal viscosity. We also measure various so-called alignment effects and discuss their influence 
on the turbulent cascade. 

Subject headings: MHD - turbulence ~ ISM: kinematics and dynamics 



1. INTRODUCTION 

Turbulence is ubiquitous in astrophysical fluids which 
are characterized by high Reynolds numbers. It af- 
fects key astrophysical processes, e.g. star formation 
(Elmegreen & Scalo 2004, McKee & Ostriker 2007). The 
observational signatures of turbulence are numerous and 
well documented. For instance, random fluctuations on 
all scales, which is a sign of turbulence, has been detected 
by a variety of observational techniques (see Crovisier 
& Dickey 1983, O'Dell & Castaneda, Armstrong et al. 
1994, Lazarian 2009). Turbulence is universal, because 
the laminar flows with high Reynolds numbers are practi- 
cal impossibility. It is driven by a variety of mechanisms 
such as supernova explosions starbursts, stellar winds, 
AGN jets, etc. In Keplerian flows turbulence is gen- 
erated by magnetorotational instability (Velikhov 1959, 
Chadrasekhar 1960, Balbus & Hawley 1991). Galactic 
disks are subject to the cosmic ray-induced Parker's in- 
stability (Parker 1966). 

Although, historically, hydrodynamic turbulence has 
been applied to astrophysics, now it is accepted that for 
almost all astrophysical fluid flows are coupled with mag- 
netic fields, at least on large scales. This necessitates 
the use of the dynamic equations that include electric 
currents and magnetic fields. The simplest approach in 
this respect is the continuous non-relativistic one-fluid 
description, known as magnetohydrodynamics or MHD. 
This approach is broadly applicable to most of the astro- 
physical environments, such as Solar wind. Interstellar 
and Intercluster Medium, molecular clouds, stars interi- 
ors, and so on, although there are some exceptions, such 
as ultra-relativistic jets and shocks, where full relativis- 
tic equations should be used or small scales of Molecular 
Clouds, where, due to relatively low ionization rate, the 
two-fluid description of ions and neutrals is more appro- 
priate. 
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The importance of astrophysical turbulence inspires 
much of theoretical and numerical work aimed at un- 
derstanding its properties. We should clarify, however, 
that there are different types of MHD turbulence. In 
this paper we deal with strong MHD turbulence, which 
can not be treated perturbatively. The theory of weak 
Alfvenic turbulence, which has limited applicability to 
astrophysics, is discussed elsewhere (Sridhar & Goldre- 
ich 1994, Ng & Bhattacharjee 1996, Lazarian & Vishniac 
1999, Galtier et al. 2000, 2002). In order to study the 
basic properties of MHD cascade and to be able to di- 
rectly compare to previous work, we restricted ourselves 
to so-called balanced MHD turbulence or turbulence with 
zero net cross-helicity. The properties of imbalanced tur- 
bulence are studied in a companion paper Beresnyak & 
Lazarian (2009). 

The issue of the spectral slopes of MHD turbulence has 
caused a substantial interest recently. A sizable num- 
ber of papers attempting to measure true asymptotic 
spectral slope of high- Reynolds number MHD turbulence 
from direct numerical simulations appeared to date. For 
example, simulations of we akly compressib l e MHD tur- 
bulence performed in Haug en et al.l (|2003l . 120041 ) used 
finite-difference code with numerical resolution of up to 
1024"^ with explicit viscous and resistive dissipation. The 
energy spectral slope for MHD was observed to be shal- 
lower than —5/3 which was interpreted as the influence 
of the bottleneck effect. A nother example is the paper 
of iMiiller fc Gr appin (2005) who measured the spectral 
slopes of decaying MHD turbulence without mean field 
and driven MHD turbulence with strong mean field. The 
pseudospectral method with ordinary viscosity was ran 
at a numerical resolution of up to 1024^. The authors ar- 
gued that the slope was close to —3/2 in the mean-field 
case. 

The motivation behind these and many other papers 
was to understand the nature of the turbulent cascade. 
The turbulent energy transfer is, in a sense, a cen- 
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tral issue of turbulence, be it hydrodynamic or MHD. 
And while hydrodynamic turbulence has its "Standard 
Model" , the nature of MHD cascade is still debated. 

An important first step in turbulence theory was made 
by Iroshnikov (1963) and Kraichnan (1965) who noticed 
that there is a local magnetic field which can not be 
excluded by a choice of reference frame, like the aver- 
age velocity in hydrodynamics. Furthermore, they as- 
sumed that turbulence is weak, because perturbations 
are smaller than the mean field. This implicitly assumed 
local isotropic dynamics and happened to be a mistake, 
since the turbulent cascade preferred perpendicular di- 
rection, i.e. produced perturbations that are more and 
more anisotropic, this way increasing interaction and pre- 
venting turbulence from becoming weak. 

The understanding of various aspects of MHD turbu- 
lence, including the role of turbulence anisotropy, com- 
pressibility, etc. resulted in a number of publications 
(see Dobrowolny, Mangeney, & Veltri 1980, Shebalin, 
Matthaeus & Montgomery 1983, Montgomery & Turner 
1984, Higdon 1984). For the most part of the paper 
we will consider incompressible MHD turbulence, which 
properties are dominated by the Alfvcnic perturbations. 
Interestingly enough, some properties of Alfvenic turbu- 
lence carry over not only to nearly incompressible low 
Mach number flows, but also for flows with Mach num- 
bers larger than unity (Cho & Lazarian 2003, see also 

It has been realized that interactions of Alfvenic modes 
in MHD turbulence has a tendency of getting stronger as 
the cascades unfolds. Goldreich & Sridhar (1995, hence- 
forth GS95) proposed a particular model of strong tur- 
bulence when the interaction strength is being controlled 
by two competing processes: a perpendicular cascade (a 
concept rigorously developed in a theory of weak Alfvenic 
turbulence, e.g. Galtier et al 2000), which tends to in- 
crease the interaction, and a decorrelation due to cas- 
cading which tends to increase the frequency of pertur- 
bations, and thus decrease the interaction. GS95 con- 
cluded that the interaction has to be marginally strong 
("critical balance"), and therefore, the cascade has to 
be of strong Kolmogorov type and has a spectral slope 
of around —5/3. Predictions of this model, such as 
scale-dependent anisotropy, were subsequently observed 
in three-dimensional simulations (Cho & Vishniac 2000, 
Maron & Goldreich 2001, etc). 

Recently, however, a number of models, motivated by 
nmncrical spectral slopes shallower than —5/3, appeared 
(Boldyrev 2005, 2006, Gogoberidze 2007). The interest 
to this field has been heated by the numerical discov- 
ery of so-called scale-dependent polarization alignment 
(Beresnyak & Lazarian 2006), which was interpreted in 
a subsequent publication by Mason et al (2006) in favor 
of Boldyrev's (2006) modification of GS95 model. 

This paper is bringing attention to serious difficulties 
that appear when one tries to measure true asymptotic 
slopes from direct three-dimensional numerical simula- 
tions which have rather moderate Reynolds numbers. 
By comparing spectra of hydrodynamic and MHD tur- 
bulence we found that MHD turbulence might be much 
more nonlocal that Kolmogorov turbulence, and, there- 
fore, require much higher resolution to obtain spectral 
slope by brute force approach. The second part of this 
paper brings polarization alignment to more numerical 



scrutiny and compares the simulation results to theory. 
Our simulations allow to test some of the existing conjec- 
tures about properties of MHD turbulence. For instance, 
these results allow us to reject a conjecture in iBoldvrevI 
12006) that the alignment is limited by the magnitude of 
local field wandering. 

In what follows in §2 we describe our approach based 
on comparing the properties of MHD and hydrodynamic 
turbulence, present our numerical setup in §3, present 
spectra in §4, discuss anisotropy in §5, discuss the in- 
teraction weakening and the alignment effects in §6 and 
§7, respectively, and provide some more hints of the non- 
locality of the MHD cascade in §8. In §9 we compare 
our numerical results with the existing theoretical pre- 
dictions, as well as with the previous numerical work; we 
also discuss the applicability of findings obtained with in- 
compressible simulations to the real-world compressible 
astrophysical turbulence. Our conclusions are summa- 
rized in §10. 

2. SLOPE MEASUREMENTS 

Various claims were made on the value of MHD spec- 
tral slope, most of which were motivated by either Kol- 
mogorov —5/3 slope of strong turbulence (Kolmogorov 
1941, GS95), or various versions of —3/2 slope (Irosh- 
nikov 1963, Kraichnan 1965, Gogoberidze 2007, Boldyrev 
2006, etc). Most numerical studies aimed to confirm ei- 
ther of the above. A number of critical issues were over- 
looked, though. Below we present a novel perspective to 
slopes measurements in MHD. It turns out that it is in- 
correct to measure slopes directly from 3D numerics due 
to a systematic error that comes with so-called bottle- 
neck effect. Also, it is impossible to measure MHD slope 
by comparison with hydro slope, because the bottleneck 
effect turns out to be different in hydro and MHD. This 
difference, however, allows us to criticize models that rely 
on strong local Kolmogorov cascading, such as GS95 or 
Boldyrev (2006) ^ 

To be precise, —5/3 is not the exact predicted slope for 
incompressible hydrodynamics. This number comes from 
the Kolmogorov self-similar cascade, but soon it was re- 
alized that realistic turbulence is not exactly self-similar. 
To correct for this, various models of intermittency were 
proposed (Kolmogorov 1962, Obukhov 1962), with the 
most popular being She-Leveque model (1994) in which 
the predicted slope is around —1.70. This slope was very 
cl ose to what was obs erved in highest-resolution DNS 
of 'Kaned a et al.l (|2003l ). In the aforementioned work it 
was possible to separate spectrum into inertial range and 
relatively more flat part that was due to a bottleneck 
effect. Although modifications of S-L model for MHD 
turbulence have been proposed, numerical studies still 
often compare the spectrum slope with —5/3. Similarly, 
when one proposes a model of turbulence with —3/2 spec- 
tral slope it is often that numerical studies aim to find 
an exact correspondence with this slope, without regard 
for intermittency. We believe that this is one important 
stumbling block in numerical determination of slopes. 

The other, probably much more important misunder- 
standing, is to disregard the systematic error that any 

^ Local cascading models normally use a formula e = pvf/ri, 
where vi is the velocity perturbation on scale r; is the cascading 
timescalo and e is a constant equal to the energy dissipation rate 
per unit volume. 
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numerical measurement of slopes in DNS brings with it. 
The bottleneck effect is a pile-up of energy before the 
dissipation scale due to the rel ative lack of en ergy in 
the dissipation range (see, e.g.. [FalkovichlfTool . Due 
to relatively low resolution of currently available simula- 
tions, this systematic error is always present. To make 
it worse, most researchers present simulations with the 
highest numerical resolution only. Although the amount 
of numerical resources available to different groups differ 
substantially, most "high resolution" simulations to-date 
have numerical boxes between 512"^ and 1024^. From 
the point of length of incrtial interval, and the influence 
of bottleneck effect, the differences in linear scale of the 
multiple of two are tiny. Another systematic error comes 
from the effect of the driving scale. Often there is a dip 
right after the driving scale, an anti-bottleneck effect of 
sorts, which appears, possibly, due to the excess of en- 
ergy on the driving scale. We are not aware of driven 
turbulence simulations that were able to get rid of this 
effect. We belive that disregarding these two effects and 
present numerical slopes as having no systematic error 
at all is wrong. 

In this paper we compare hydrodynamic and MHD en- 
ergy slopes obtained with the same code, the same driv- 
ing and exactly the same linear dissipation. Since there 
are good theoretical predictions for asymptotic isotropic 
hydro turbulence, we can try to use those. If one finds 
that the nature of energy transfer in MHD and hydro is 
similar (this is suggested in GS95 model where a strong 
local Kolmogorov-like cascade is assumed), then we can 
directly compare MHD and hydro slopes and make state- 
ments on MHD slope. Unfortunately, as we show in 
the two subsequent sections, this is not the case. The 
defining feature of our simulations is the use of different 
types of linear dissipation, namely natural viscosity and 
hyper-viscosity. Although there had been some similarity 
in spectral slopes of MHD and hydro in normal-viscous 
case, the hyper-viscous cases were very different. This 
suggests that the nature of MHD and hydro cascades are 
different and one can not use slope comparison between 
MHD and hydro to get rid of the aforementioned system- 
atic error. 

3. NUMERICAL SETUP 

Incompressible MHD and Navier-Stokes equations can 
be written in the following simple form 

dtw^ + 5(wT • V)w± = -M„(-V2)"w±, (1) 

where S* is a solenoidal projection and w='= (Elsasser 
variables) are defined in terms of velocity v and magnetic 
field in velocity units b = B/(47r/9)^/^ as w+ = v -f b 
and w~ = V — b. Navier-Stokes equation is a special 
case of equations (1), where 6 = and, therefore, both 
equations are equivalent with = w~ . The RHS of 
this equation is a linear dissipation term which is called 
viscosity or diffusivity for n = 1 and hyper-viscosity or 
hyper-diffusivity for n > 1. Here we assumed that vis- 
cosity and magnetic diffusivity is the same for velocity 
and magnetic field. This is almost never true for re- 
alistic astrophysical plasmas. However, as long as one 
wants to study the dynamics of large scale turbulence, 
it is acceptable. This is because the dissipation terms 
in today's numerical simulations are never as small as to 
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TABLE 1 

Description of simulations. * - At is the duration of the 

HIGH RESOLUTION RUNS, PRIOR TO THAT WE RAN LOWER 
RESOLUTION RUNS FOR A LONG TIME TO ENSURE STATIONARY 
STATE. E.G. WE RAN Bq = CASE FOR 180 TIME UNITS PRIOR TO 
GOING TO HIGH RESOLUTION. ONLY LAST 3-4 TIME UNITS OF HIGH 

RESOLUTION RUNS WERE USED FOR MEASUREMENTS. ONE TIME 
UNIT CORRESPONDS TO AN EDDY TURNOVER TIME OF THE LARGEST 
COHERENT EDDIES. 

simulate real physical dissipation, instead, they are used 
to remove energy on small scales and assure stability of 
the code. In finite-difference codes a numerical dissipa- 
tion is always present and the linear dissipation terms 
are often omitted altogether (see, e.g.. Stone et al 1998). 
In pseudospectral code, such as our own, the energy is 
conserved with rather good precision, so the use of linear 
dissipation is necessary. 

We evolved incompressible MHD and Navier-Stokes 
equations in time using a well-known pseudospectral 
technique (see, e.g.. lCho fc Vishnia^l200CI[ ). We have cho- 
sen pseudospectral code as it allows precise control over 
dissipation. Our hydro code was identical to the MHD 
one, except, naturally, for the lack of magnetic field. The 
summary of high-resolution runs is presented in Table 1. 

We performed four types of simulations: a simulation 
of statistically isotropic hydro turbulence; a simulation 
of well-developed stationary MHD turbulence without 
mean field {Bq = 0), which also has been called statisti- 
cally isotropic MHD turbulence; a simulation of so-called 
transAlfvenic turbulence (with Bq = 1, SB Sv ~ 1); a 
simulation of strong anisotropic subAlfvenic turbulence 
(with Bq = 10, 6B Sv ^ 1). Special attention had 
to be taken to the last case, where, in order for the 
turbulence to be strong, the fields had to be strongly 
anisotropic on the outer scale. To ensure this, the nu- 
merical box was elongated in real space and x-direction 
(a direction of the strong mean field) was 10 times longer 
than y and z directions. The driving had anisotropy that 
corresponded to the dimensions of the bo x. Thus in this 
case our setup is similar to one in Maron fc GoldreichI 
(j2001). 

We used random solenoidal velocity driving in k-space 
between k=2 and 3.5. The largest coherent eddy size L 
(determined by structure function technique) was around 
1/4 of a box size (27r) and the eddy turnover time for 
this scale was around unity. The Alfvenic self-crossing 
time for this scale was also around unity ^. The self- 
correlation timescale for our driving force was t = 2 for 
all wavemodes. This is somewhat longer than the eddy 
turnover time. We performed a test study of the force 
correlation time influence on the self-correlation time of 

2 In the strong field case, Bq = 10, the Alfven speed was ten 
times higher, but the largest coherent eddy was elongated with 
x/y aspect ratio of around 10, so this gives the same estimate for 
Alfvenic crossing time. 
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velocity, taking r = 1/4, 2, 8 and found no systematic de- 
pendence. The self-correlation timescale of velocity was 
always around unity. We concluded that, in the case 
of strong turbulence, the self-correlation timescale of ve- 
locity is primarily determined by nonlinear interaction, 
rather than driving . We ran another, higher resolution 
simulation, trying to find a difference in slopes between 
simulations with r = 2 and r = 6 and found none. Note, 
that n umerical studies do vary in terms of t. For in- 
stance, iMiiller & Graooin ('2005") used driving constant 
in time, and Haugcn ct al. (2004) used (5-correlated driv- 
ing. Our tentative conclusion is that the difference in r 
marginally affects the results of numerical simulations of 
turbulence. In addition, in simulation M5 we used El- 
sasser driving ^ instead of velocity driving, to check if it 
changes the tendencies observed in previous runs. 

One of the advantages of driven versus decaying simu- 
lations is that it describes a stationary random process, 
so, by applying ergodic hypothesis one can approximate 
statistical averages with time-averages. While, in decay- 
ing simulations, if one wants to average over time to re- 
duce fluctuations, some hypothesis on the decaying pro- 
cess has to be adopted ^. 

4. SPECTRA (3D- AVERAGED AND ID). 

Usually, bottleneck effect, a pile-up of energy near 
the dissip ation scale, is di scussed in relation to hyper- 
viscositv (ICho fc Vishnia^ [2000) or numerical viscosity 
(|Kritsuk et al.l 2008f). But t he bottleneck effect was 
also predicted (lFalkovichlll994l ) and numerically observed 
(jKaneda et al.l l2003[ ) for normal viscosity. This means 
that the slopes, measured in limited resolution simula- 
tions with normal viscosity, do not necessarily reflect true 
asymptotic slopes. 

Fig. 1 shows Ek , a three-dimensional power spectrum 
F(k.) integrated over angle in /c-space: 



F(k)dk. 



(2) 



|k|=fc 



This is the most popular choice for spectrum in numer- 
ical turbulence since it is fairly straightforward to cal- 
culate. Because fc-space is discrete, the integration is a 
summation of energies of all modes that lay in a shell of k- 
magnitudes between n and n+1. The number of modes, 
that fall into each shell, fluctuates, so this spectrum has 
a characteristic toothed shape which is not remedied by 
time integration. We plot Ej. so that our results can be 
directly compared to previous numerical work that used 
this quantity. 

In the strongly anisotropic, subAlfvenic simulations (as 
in runs M4 M5) it makes sense to measure a perpendic- 
ular spectrum (i.e., integrated over and then over the 
direction of k±). This spectrum, however, was virtually 
identical to 3D spectrum integrated over solid angle, i.e., 
Ek, so, for uniformity, we plot only Ek in all cases. 

^ One can argue that in case of the absence of nonlinear term 
the equation dtv = f will give infinite correlation time for v even 
if / self-correlation time is finite. 

* Each of the pair of Elsasser-field equations in MHD bear close 
resemblance to Navier-Stokes equation. Elsasser variables are de- 
fined as = V ± b, where b is magnetic field in velocity units. 

^ Normally, one wa nts to normalize the spectrum as in 
IMiiller fc GrappinI II2005I) . but this entails the hypothesis that spec- 
tra are similar at different stages of decay. 
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Fig. 1. — Three-dimensional angle-integrated spectra for all runs, 
compensated by a factor of k^^'^. We multiplied the spectra by a 
quotient of 3 to separate them from each other. Top: HI, Ml; 
middle: H2, M2, M3, M4; bottom: H3, M5. A significant increase 
in perceived "inertial interval" between top plot and bottom plots 
is mostly due to difference between normal viscosity and hyper- 
viscosity, rather than resolution. 

The parallel spectrum (integrated over k±) was of no 
interest to us, since there are no clear theoretical predic- 
tions for it ^. 

We also measured so-called one-dimensional spectra 
Pk , which is a power spectrum of the vector field sampled 
along an arbitrary line and then averaged over ensemble. 
It can also be written as a Fourier transform of the cor- 
relation function B{r) (Monin & Yagiom 1975): 



-ikr 



B{r)dr. 



(3) 



Although most simulations measure iJj., the struc- 
ture/correlation function scalings are the primary predic- 
tions of the Kolmogorov model, so it makes more sense to 
measure Pk rather t han Ek- Also. Pk is less prone to bot- 
tleneck effect (Do bler et al.l [20031 ). We measured Pk by 

^ One can relate this spectrum to the parallel structure func- 
tion, calculated with respect to global mean field. This SF, how- 
ever, do not show properties, predicted by GS95 model. The pre- 
ferred way is to calc ulate SFs with respect to local mean field 
llCho fc VishniadlWa) . 
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Fig. 2. — One-dimensional spectra, compensated by fc^/^. 

averaging over directions in each particular snapshot and 
then averaging over time. Pfc is presented in Fig 2. For 
the fully statistical homogeneous isotropic sample with 
infinite spacial resolution and infinit e dimensions one 
can d erive relation — —kdPk/dk (jMonin fc YaglomI 
|1975( ). Numerically speaking, for a discrete sample with 
periodic boundaries this relation is satisfied fairly well. 
It ensures that in an infinite inertial interval Ek and Pk 
will have the same slope. In a limited resolution of our 
simulations the slopes of Ek and Pk are fairly different 
though. We feel that this difference is another useful in- 
dication that numerical determination of slopes is rather 
limited (see also §2). 

5. ANISOTROPY 

Hydrodynamic turbulence is presumed to be isotropic 
on small scales, while MHD turbulence is anisotropic and 
this anisotropy increases to small scales without limit 
(GS95). The scale-dependent anisotro py, pr e dicted by 
GS95 was first observed in iCho fc Vishmag (pOOOl ) by 
the method of second-order structure functions and con- 
firmed in a number of subsequent publications. It is crit- 
ical that structure function is calculated with respect to 
the local magnetic field, otherwise the scale-dependent 
anisotropy is not observed (see discussion on various def- 
inition of the "local mean field" in Beresnyak fc Lazarian 
2008). Fig. 3 shows two-dimensional structure function 
where x-axis measured distance along local mean mag- 
netic field while y-axis measured distance perpendicular 
to the field. In hydrodynamic case the "preferred" direc- 
tion was chosen arbitrarily. Quite predictably, the hydro- 
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Fig. 3. — Contours of the structure function < |ti(r-(-i?) — ?;(r)|^ -|- 
\b(r + R) — f'(r)p >r where R takes arbitrary angles with respect 
to the local magnetic field (for hydro case this direction was chosen 
to be along x). X-axis is perpendicular component, while y-axis is 

2/3 

parallel component. Dashed line correspond to GS95 ry ~ . 
Hydro turbulence show isotropy on small scales, while MHD tur- 
bulence show scale-dependent anisotropy. 

dynamic turbulence is isotropic. The MHD turbulence, 
however, is scale-dependently anisotropic. For example, 
Bq — Q case is isotropic on outer scale and Bq ~ 1 case 
is almost isotropic on outer scale, but on small scales 
both are anisotropic with anisotropy ratio of around 4. 
The Bq = 10 case is anisotropic on outer scale and this 
anisotropy increases by about a factor of 4 towards small 
scales. These results are approximately consistent with 
the predictions of GS95. 
In order to study deviations from prediction of GS95, 

ry ~ r]l^^ , we created a correspondence between parallel 
and perpendicular scales and plotted it on Fig. 4. This 
correspondence was achieved by finding equal values of 
parallel ar id perpendicular second o rder structure func- 
tions a s in lCho fc Vishniad (|2000D or lMaron fc GoldreichI 
(|2001h . The observed deviations could be connected to 
nonlocality, and/or alignment effects which are discussed 
further. 

6. NUMERICAL EVIDENCE OF INTERACTION 
WEAKENING 

In GS95 turbulent non-linear interactions are strong 
and the cascading happens fast. In this approach the 
Kolmogorov's ^ —5/3 spectral slope is expected. In or- 
der to explain shallower slopes Boldyrev (2006) conjec- 
tured that interaction is depleted in strong anisotropic 
turbulence by a factor which is similar to that in the 
Iroshnikov-Kraichnan (IK) model, i.e. Sv/va- In sim- 
ulations with real viscosity we obtain MHD dissipation 
scale d^MHD which is somewhat larger than hydro dis- 
sipation scale di/HD, which can be an indication of de- 
pletion of interaction. Figs 1 and 2, upper panel, indi- 
cate that di/MHD is approximately 1.3 times larger than 
c^i^HD, which is consistent with the difference in dissipa- 
tion scales of Kolmogorov and IK models, assuming the 
length of the inertial interval L/d^ of around 7. However, 
this result can also be explained by the difference in Kol- 
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Fig. 4. — Deviation of anisotropy from GS 
law IGoldreich fc jjridhaj JTMST) A ~ A^/S. M2, M3 and M4 
are shown. M5 shows behavior similar to M4, and Ml is similar 
to M3. 

mogorov constant Ck of hy dro and MHD turbulence. 
Indeed, for MHD Ck ~ 2.2 (jBiskampI [200l . while for 

hydro Ck ^1-6 (|Frischlll995f ). also note that ~ C^J^. 
Thus, this way of proving the weakening of interaction is 
still controversial. 

7. ALIGNMENT EFFECTS. 

While most MHD turbulence models use mean-field 
approach and assume that turbulence is characterized 
fully by spectrum or structure functions of the fields, 
i.e. {6v'^) and {Sb'^), recently a considerable attention has 
been drawn to so-called alignment effects (iB oldyrcv 2005; 
iBeresnvak fc Lazarianll2006l : lBoldvrevll2006D . The align- 
ment effects can be understood, in general, as a property 
of multi-variate pdf of the fields containing various corre- 
lations. The scale- independent alignment effects are not 
so interesting, because they can only modify the Kol- 
mogorov constant of turbulence, while scale-dependent 
alignment can, in principle, modify the slope. 

Consider the alignment of Alfven mode when all per- 
turbations are perpendicular to the local magnetic field, 
i.e. lie in the same plane. For this purpose we use struc- 
ture functions where vectors are projected on 1 x B di- 
rection where 1 is the direct ion, connecting two points 
(|Beresnvak fc Lazarian|[2"006D . While a study a full mul- 
tivariate PDF could be an overwhelming task, one can 
introduce a few statistical measures of alignment that 
could be of interest. Alignment, derived from zeroth or- 
der statistical moments: angle alignment, 

AA= (|sin6'|), (4) 

where 9 is an angle between Elsasser variables perturba- 
tions (5w+ = (5v + (5b and Sw^ = (5v — (5b, another angle 
alignment, 

AA2 = (|sin^2|), (5) 

where 02 is an angle between (5v and Sh. Alignment 
derived from second o rder statistical mo ments, "polar- 
ization intermittency" (jBeresnvak fc Laz arian 2006): 

PI = {\6w+6w- sin6'|)/(|(5w+(5ii;"|); (6) 



imbalance measure 

IM = {\diw+f - diw-f\)/(6{w+f + 6iw-f); (7) 

imbalance correlation 

IC = {\5w+Sw-\)/{6{w+f + 5{w-)^), (8) 

velocity and magnetic fields correlation 

FC ^ {\Sv\\5b\)/{5v^ + Sb^) (9) 

(we found that FC is almost constant in our measure- 
ments and it will be assumed constant thereafter). IM 
and IC are describing the same effect, namely dynamic 
imbalance between Elsasser variables dw~^ and S'w~ , but 
in a different way. For independently distributed Gaus- 
sian fluctuations one has AA = AA2 = PI = IM = 2/t:, 
IC = I/tt. 

Physically, one may have only two types of alignment, 
because there are four variables (two vectors) minus 
normalization and minus one arbitrary rotation along 
axis perpendicular to the magnetic field. Let us choose 
PI and IC as two measures of alignment. PI is inter- 
esting, because it is the factor, by which an interac- 
tion is reduced in a nonlinear MHD term {Sw ■ W)Sz ~ 
{Sw • kz)dz ~ Swkz Sz sin with re spect to mean field es- 
timate of SwkzSz ()Boldvrevll2005[ ). This quantity, along 
with AA, w as first measured in numerical simulations by 
IBeresn v ak fc Lazarian (2006 ). In the subsequent publi- 
cation lMason et al.l (|2006h used a second order structure 
function measure, very similar to PI, termed "dynamic 
alignment" 

DA= (|,5u(5&sin6i2|)/(|(5w||(56|). (10) 

It can be expressed, to a constant, as DA ^ PLIC/FC ^ 
PI • IC (since w+ x w~ = — 2v x b), i.e. it is a combi- 
nation of polarization intermittency and imbalance cor- 
relation IC. It is not clear yet, whether intermittent im- 
balance can reduce interaction in the balanced turbu- 
lence. For example, if one es timates energy dissip ation 
as Sw+Sw~{Sw+ + Sw~)/X (|Lithwick et alJ 120071 ). it is 
insensitive to the dynamic imbalance 6w^ ± e, to the first 
order of e. But one also may argue that the imbalanced 
case is more complicated and the i nteraction is reduced 
(IBeresnvak fc Lazarianll2008l . l2009f) 

Numerically speaking, alignment effects were very sim- 
ilar for all sub-Alfvenic and trans- Alfvenic cases'^. Fig. 
4 shows different alignment measures described in this 
section for M5 {Bq = 10, Elsasser driving). In the mid- 
dle of the inertial interval alignment factors depend on k 
as AA ~ A:0 022, aA2 ~ fco osa, pi ^ /fcO.ios^ ~ k°-°^\ 
IM - DA - A:0-207. Note, that DA is short 

of the fcO-25 dep endence predicted for "alignment an- 
gle" in iBoldvrevI ([2006). The alignment dependence in 
the velocity-driven, B = 10 M4 simulation is somewhat 
different: AA ~ fc0.028^ aA2 ~ fc0.052^ pj _ ^o.i37^ 

IC - fc°■05^ IM ~ DA ~ If one wanted 

to explain the difference in slopes between M4 and M5 
(~ 0. 05. Fig 2) by the difference in DA slopes (accord- 
ing to IBoldvrevI (j2006[ ) the DA slope a adds 2/3a to the 
spectral slope), such an explanation would be impossible. 

^ Beresnyak & Lazarian (2006) observed alignment not only in 
incompressible simulations, but in trans- Alfvenic, trans-sonic com- 
pressible simulations as well. 
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Fig. 5. — Various alignment measurement from M5. 



On the other hand, M4 and M5 are significantly differ by 
IM, which suggests that shallow slopes are primarily due 
to imbalance. 

M3 {Bq = 1, velocity driven) simulation show align- 
ment effects which are simila r to M4. This di rectly con- 
tradicts to the prediction in 'B oldvrevI (j2006f ). that the 
alignment is a function of 6B/Bq, as we observe essen- 
tially the same alignment in simulations with SB^/Bq = 
1 and SBl/Bq — 0.1 (see also Fig. 6). We therefore con- 
clude that the alignment is an intermittency effect^ that 
accumulates along the cascade, rat her than bein g deter- 
mined rigidly by SB\/Bq as in Bol dvrevI (|200l) . Alter- 
natively, the alignment could be a natural feature of the 
nonlocal energy transfer. This calls for further investiga- 
tion. 

8. OTHER HINTS ON NONLOCALITY. 

In the sub-Alfvenic turbulence, the properties are well- 
represented by Elsasser variables and, if we assume 
locality of interaction of Sw^ and 5wf , the properties 
of 6vi and 6bi are supposed to be identical. However, 
we observed notable differences even after large statisti- 
cal averaging. Namely, the magnetic energy was higher 
than kinetic energy in subAlfvenic runs M4 and M5. 
The last case, which was driven by Elsasser variables, 
show 50% magnetic energy excess on large scales and 
around 20% on small scales. This so-called residual en- 
ergy which was proposed to have "-2" spectrum scaling 
(IMiiller fc Grap pin 2005|) but was somewhat shallower 
than "-2" in this Elsasser driven run. Other runs did 
not show any regular scaling for residual energy. For ex- 
ample, statistically isotropic MHD turbulence M2 show 
dominance of kinetic energy on large scales (which is typ- 
ical for MHD turbulence driven with velocity on outer 
scale), but on small scales magnetic energy dominates 
(which is typical for almost any driven MHD turbulence) . 
This reinforces our conjecture that currently available 3D 
MHD simulations do not yet exhibit inertial ranges and 
the flat portion of the spectrum can not be considered a 
part of the inertial range of local Kolmogorov-type tur- 

® An intermittency eorrection is usually written as a function of 
the ration of the scale in question to the outer scale, e.g. Svf = 

C (el)'^ ^ ^ (I / L)" , where {l/L)°' is the intermittency correction. In 
our case we also see that the alignment is a function of l/L. 



bulence until a number of other conditions are satisfied, 
among which is an equipartition between spectral kinetic 
and magnetic energies. 

9. DISCUSSION 
9.1. Theory 

The nature of the turbulent cascade and the slope of 
the spectrum of fluctuations is a central issue of tur- 
bulence and has been discussed in a majority of pa- 
pers devoted to turbulence theory and numerics. The 
hydrodynamic isotropic turbulence has its "Standard 
Model" which is based on Kolmogorov's assumption of 
self-similarity ^ and Kolmogorov's "-4/5 law": 

{{6v\i{l)f) = -^d. (11) 

Assuming similarity {{Sv\\{l))'^) ^ {{6v{l)y}'^/^ one ob- 
tains the second order structure function slope of 2/3 
which correspond to spectral slope of —5/3. In MHD, re- 
lations, similar to Kolmogorov's "-4/5 law", exists, e.g., 

{Swjil)i6w^{l)r)^-^e^l (12) 

(Chandrasekhar 1951, Politano & Pouquet 1998)^°. 
However, in MHD case, this does not directly hints on 
scalings for {{Sw^ (l))'^) and it is not clear which similar- 
ity hypothesis has to be adopted. A nice demonstration 
of this is to apply the aforementioned exact relations to 
the case of weaA: turbulence, where, in the ansatz of three- 
wave interaction, one has to obtain {SwJ {l){Sw^ {l))"^) ^ 

a{l){{6w{l))'^)/vA (anisotropy a{l) is defined as a ratio 
of parallel scale A(/) to perpendicular scale I), which 
is very much unlike the Kolmogorov similarity hypothe- 
sis. Thus, the spectral slope scalings are still uncertain, 
which stimulates further research in this field, including 
attempts to measure the slope numerically. 

The nonlocality of turbulent energy transfer has been 
claimed in quite a few publications. In the discussion 
of this numerical work we will mention the most rele- 
vant ones and defer more exhaustive discussion to future 
review papers. 

It was suggested in iGogoberidzd ()2007| ) that in MHD 
turbulence with strong mean field the nonlinear inter- 
action could be nonlocal, with outer scale perturbations 
decorrelati ng high-frequency interactiiig edd ies, leading 
to IK-type (jIroshnikovl(l963l : lKraichnanll965l ) interaction 
weakening and —3/2 spectrum. However, in aforemen- 
tioned model the energy cascading itself is local and 
is performed by high-frequency eddies, i.e. there is no 
energy transfer from large scales directly to small scales. 
We conclude that this model can not explain the lack of 
bottleneck effect, observed in simulations. 

In a recent model o f imbalanced MHD turbulence 
(jBeresnvak fc Lazarianlf2008i ) the eddies of the dominant 

^ This assumption is not precisely satisfied because of so-called 
intermittency. The most popular descriptions of intermittency are 
probabilistic models (Obukhov 1962, Kolmogorov 1962, She & Lev- 
eque 1994). They give corrections to spectral slopes and higher- 
order scaling exponents assumed in self-similar model 

D = 3 for statistical isotropy and isotropic structure function 
or D = 2 for cylindrical symmetry and perpendicular structure 
function, e.g. in turbulence with strong mean field. 
In a sense that e = /ri is used, see also §2. 
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Fig. 6. — A comparison of DA in sub-Alfvenic (M4, dotted) and 
trans-Alfvenic (MHDSblh, dashed) simulations with &B/B pre- 
diction of iBoldvrevi (i2006i ). &B/B is estimated from 2nd order 
transverse (or isotropic) structure function. 

component on a certain scale are aligned not with re- 
spect to the loca l mag netic field on the same scale as in 
iCho fc Vishniad ()2000D (balanced turbulence), but with 
respect to magnetic field on some larger scale. One 
may speculate that even in the balanced case the dy- 
namic imbalance can cause polarization alignment, or, 
more likely, that thes e eff'ec ts are interrelated. While 
iBeresnvak fc LazarianI ()2008D . by itself, is a mean field 
model and reproduce locality, "-5/3" spectrum and "2/3" 
anisotropy of GS95 in the balanced limit, its future exten- 
tions to include local imbalance and polarization align- 
ment seems promising. 

Does the hypothesis of Boldyrev (2006) that the v and 
b alignment is limited by field wandering is justified from 
theory ground? In a little thought experiment one can 
imagine a perfectly aligned state where the magnitudes 
of V and b are equal. This is a case of a perfect im- 
balance and also an exact solution of MHD equations. 
Such a state will propagate without distortion, in other 
words, no de-alignment is going to happen, although this 
state certainly has some level of field wandering. This 
thought experiment hints that the hypothesis of v and 
b alignment being limited by field wandering directly 
contradicts MHD equations. The effect of local imbal- 
ance, therefore, has to be treated in a more complicated 
way. Fig. [B]shows a comparison between field wandering 
{5B/B) and "dynamic alignment" in transAlfvenic M3 
and subAlfvenic M4. We see that in subAlfvenic case 
the dynamic alignment is an ord er of magnitude off the 
SB/B, which is in contrast with IBoldvrevi (|2006( ) which 
suggest that DA 5B/B. 

9.2. Previous numerical work 

The study o f weakl y coru pressible MHD turbulence by 
iHaugen et gdl (|2003l . (2004( 1 revealed some difference in 
bottleneck effect in hydro and MHD cases, although the 
authors reported bottleneck behaviour as "very similar" . 
The difference was small, possibly, due to adopted first 
order viscous and resistive dissipation. Aforementioned 
work, however, used finite-difference code and, inadver- 
tently, the numerical dissipation had a different character 
in MHD and hydro simulations, which precluded a rig- 
orous comparison. Nevertheless, we can say that it is 
consistent with what we see in our incompressible simu- 
lations. 

Biskamp et al. (1998) studied two-dimensional MHD 



and EMHD turbulence. They noticed that both MHD 
and EMHD cases have an unusual nonlocal bottleneck 
effect, which appeared differently depending on numer- 
ical resolution. It is possible that spectral fiattening 
from such an effect can be perceived as a "false" in- 
ertial interval and lead to an incorrect estimate of the 
slope. Although there had been much discussion on the 
analogy between two-dimensional MHD turbulence and 
MHD turbulence with strong mean field, the predom- 
inant picture is that Alfvenic turbulence is essentially 
three-dimensional (GS95, Cho & Vishniac 2000, Maron 
& Goldreich 2001, Cho et al. 2002). 

Yousef et al. (2007) claimed that in the statistically 
isotropic MHD turbulence (similar to our simulation M2) 
one has a folded magnetic field structures that directly 
non- locally interact with outer-scale motions. They con- 
cluded that one can not use Kolmogorov argumentation 
because of this nonlocality. We, however, belive, that 
Kraichnan's argumentation regarding the dominance of 
local mean field, i.e. the Alfven effect, is correct even in 
the case of B^ = 0. This is somewhat hard to demon- 
strate in 3D numerics, however, because for Bq = Q there 
is a transition region between outer scale and inertial in- 
terval where kinetic energy dominates and the magnetic 
spectrum is very fiat, i.e., there is no clear dominance of 
the large-scale magnetic field. Also the use of first-order 
(natural) viscosity in aforementioned three-dimensional 
simulations made the inertial range very short, which 
created an illusion of a universal folded magnetic field 
structure. 

Alexakis et al. (2005a, 2005b) used a specific numeri- 
cal tool to quantify the transfer of energy between scales 
in both hydro and MHD turbulence. They claimed that 
MHD case is somewhat nonlocal. However, the more 
radical claim is that the hydro cascade and the larger 
part of MHD cascade is extremely local, i.e. the energy 
is not transferred between k and 2k, as in Kolmogorov 
model, but instead between k and k + ko where ko is de- 
termined by outer scale, which breaks self-similarity. We 
find these results rather puzzling and defer discussion un- 
til independent confirmation is available. Until then, we 
assume that Kolmogorov picture is roughly appropriate 
for isotropic hydrodynamic turbulence with large inertial 
range. 

9.3. Implications for astrophysical turbulence 

Realistic astrophysical turbulence is, in general, com- 
pressible. The examples of weakly compressible flows 
are the quiet convection in main sequence stars and tur- 
bulence in very hot intracluster gas. The turbulence in 
most of the Interstellar Medium, however, is strongly 
compressible, due to effective cooling. This raises the 
question of to what extend the results of incompressible 
simulations are applicable to astrophysics. Intuitively, 
one could expect that weak small-amplitude Alfven and 
slow-mode turbulent perturbations should be nearly in- 
compressible. The question is how to deal with the fast 
mode and large-amplitude MHD turbulence in compress- 
ible fluids. 

The coupling of Alfvenic motions and compressible mo- 
tions is a difficult subject. Theoretical arguments that 
fast, slow and Alfven modes may create independent en- 
ergy cascades were provided in GS95, Lithwick & Goldre- 
ich (2001) and Cho & Lazarian (2003). These arguments 
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were applicable to strong Alfvenic turbulence^^. Cho & 
Lazarian (2002, 2003) numerically demonstrated that for 
even for appreciable Mach numbers (up to 10), the prop- 
erties of Alfvenic modes (scalings and anisotropy) were 
similar to those in incompressible simulations^^. They 
also showed that the decay time for the Alfvenic modes 
is fast and not mediated through coupling with com- 
pressible motions, which used to be the common wisdom 
at the time of the study. The effect of scale-dependent 
polarization alignment, a characteristic of Alfvenic cas- 
cade discussed in this paper, happen to be present in 
both compressible and incompressible MHD simulations 
(Beresnyak & Lazarian 2006). While the extend to which 
strong shocks can modify the Alfvenic cascade deserves 
more study, the arguments above make us confident that 
the studies of incompressible turbulence are of primary 
importance to understand astrophysical turbulence. 

There is another issue, which is frequently ignored 
when one compares numerical MHD with interstellar 
turbulence. If magnetic fields are perfectly frozen into 
the fluid, turbulence creates multiple small-scale current 
sheets which arc difficult to dissipate. Within such zones 
the frozen-in condition is no longer valid and magnetic 
reconnection takes place (see Biskamp 2000, Priest & 
Forbes 2000, Bhattacharjee 2004, Zweibel & Yamada 
2009). The Lundquist number S = LVa/ti that char- 
acterizes how well magnetic fields are frozen in (L is the 
scale of the current sheet and rj is magnetic diffusivity) , 
is very high, e.g. > 10^°, for most astrophysical flu- 
ids and is fairly low, e.g. < lO"*^, for MHD simulations. 
If magnetic reconnection in astrophysics depends on S, 
this presents not only a problem for most of MHD sim- 
ulations, e.g. simulations of dynamo, molecular clouds, 
accretion disks, but also means that the numerical re- 
sults on turbulent scalings may not be trusted. We be- 
lieve that the extensive observational data suggests that 
magnetic reconnection is fast. Also, a model predicting 
fast reconnection in turbulent fluids, Lazarian & Vish- 
niac (1999), has been successfully tested in Kowal et al. 

The more detailed study (Chandran 2005) of weak Alfvenic 
turbulence in a pressure-less medium showed that the interactions 
between the fast and Alfvenic modes may be significant in a certain 
regions of k-space. 

Density perturbations, which are absent in incompressible 
flows, are definitely affected by compressible motions. However, 
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(2009). This gives additional support to numerical test- 
ing of astrophysical turbulence. 

The slopes and turbulence anisotropics are important 
for a variety of astrophysical phenomena. For instance, 
scattering and turbulent acceleration of cosmic rays de- 
pends on the scaling of MHD turbulence (see Chandran 
2000, Yan & Lazarian 2002, 2004). So does the per- 
pendicular diffusion of cosmic rays and heat transport 
in plasmas (Narayan & Medvedev 2001, Lazarian 2006). 
Naturally, it is important to establish the true scalings 
of MHD turbulence. This paper testifies that a higher 
resolution numerical simulations are required for accu- 
rate testing of the present and future MHD turbulence 
models. 

10. CONCLUSIONS 

We conclude that although the Kolmogorov-like 
Goldreich-Sridhar (GS95) model is appealing, simple and 
captures some essential physical properties of the strong 
MHD turbulence, such as scale-dependent anisotropy, it 
should be amended to explain cascade nonlocality and 
scale-dependent alignment effects. How to achieve this 
is the issue of the future research, as we demonstrate that 
the existing attempts to improve GS95 do not agree well 
with the presented numerical simulations. In addition, 
we issue a note of warning that the numerical measure- 
ments of the spectral slope that served as a motivation for 
many of theoretical studies are unlikely to represent the 
true theoretical slopes due to the non-locality of MHD 
cascade. 
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Beresnyak et al. (2005) found that the structure of logarithm 
of density reveal scale-dependent anisotropy, similar to GS95 law. 
This presumes that a significant fraction of density structures are 
created by shearing by Alfvenic perturbations. 
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